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We present a new method, Non-Stationary Forward Flux Sampling, that allows efficient simulation of rare 
events in both stationary and non-stationary stochastic systems. The method uses stochastic branching and 
pruning to achieve uniform sampling of trajectories in phase space and time, leading to accurate estimates for 
time-dependent switching propensities and time-dependent phase space probability densities. The method is 
suitable for equilibrium or non-equilibrium systems, in or out of stationary state, including non-Markovian or 
externally driven systems. We demonstrate the validity of the technique by applying it to a one-dimensional 
barrier crossing problem that can be solved exactly, and show its usefulness by applying it to the time- 
dependent switching of a genetic toggle switch. 



I. INTRODUCTION 

Rare events, which occur infrequently but have im- 
portant consequences, control the dynamical behavior of 
many physical systems, both in and out of equilibrium - 
classic examples include crystal nucleation, protein fold- 
ing, earthquakes and traffic jams. When simulating such 
systems on a computer, some form of enhanced sam- 
pling is usually needed in order to generate any signif- 
icant number of rare event samples on the time scale of 
the simulation. While a number of enhanced sampling 
methods are available for systems in steady state, many 
important rare event processes happen in non-stationary 
systems, for which most existing methods are unsuitable. 
In this article, we introduce a new enhanced sampling 
scheme, Non-Stationary Forward Flux Sampling, which 
allows efficient simulation of rare events in both station- 
ary and non-stationary stochastic systems. 

For systems in thermodynamic equilibrium, a large va- 
riety of rare-event techniques h ave been developed. One 
is the Bennett-Chandler method^, which involves a cal- 
culation of the free energy along a predetermined reac- 
tion coordinate, followed by a computation of the ki- 
netic pre-factor by firing off trajectories from the top of 
the free-energy barrier. Other techniques are transition 
path sampling^ and transition interface sampling (TIS) 4 , 
which employ Monte Carlo sampling of the ensemble of 
transition paths, approximate schemes such as partial- 
path TI9^ and milestoning 6 -!, which use a series of inter- 
faces in phase space between the initial and final state, 
and string methods-^. While these schemes have been 
successfully applied to a large class of problems, they do 
require knowledge of the phase-space density, which lim- 
its their use to systems in thermodynamic equilibrium. 

For non-equilibrium systems, the phase-space density 
is generally not known. This severely limits the possibili- 
ties for devising enhanced sampling schemes to calculate 
transition rates. Yet, for non-equilibrium systems that 
are in stationary state, recently a number of rare-event 
techniques have been developed. One is the weighted- 
ensemble (WE) methocP^, where phase space is divided 
into bins, and trajectories are selected and re-weighted 



bin- wise to achieve uniform coverage of the phase space. 
Another technique is the non-equilibrium minimum ac- 
tion method^, which allows the characterization of tran- 
sition paths but not rate constants. Non-equilibrium um- 
brella sampling^ coarse-grains systems with Markovian 
dynamics on overlapping grids in state space and biases 
inter- vs. intra-bin transitions. Forward-Flux Sampling 
(FFS^P^uses a series of interfaces in phase space between 
the initial and final state to drive the system over the 
barrier in a ratchet-like manner, by capitalizing on those 
fluctuations that move the system from one interface to 
the next. While these methods do not require thermal 
equilibrium, they rely on the system being in stationary 
state. 

In reality, however, many important rare event pro- 
cesses happen in systems which are not in stationary 
state. For these processes, the propensity (probabil- 
ity per unit time) for the rare event to occur is time- 
dependent; this time dependence may be caused by ex- 
ternal driving, by transient relaxation of the system from 
an out-of-equilibrium initial state, or by the presence of 
memory in the dynamics on a relevant time scale. In 
fact many real-life instances of the rare event processes 
mentioned above are time-dependent, such as: crys- 
tal nucleation during flash-freezing (e.g. when preparing 
cryo-electron microscopy samples); protein folding dur- 
ing transient association with a chaperone protein; and 
triggering of traffic jams by brief disturbances on the 
road. Other interesting cases include transitions between 
multiple limit cycles in neural networks under time- 
dependent stimuli (as suggested for epileptic seizures, 
e.gP^), and the response of metastable biochemical net- 
works to transient signals, e.g. in cell differentiation^ or 
in viral life cycle progression, see sec. |IV| These impor- 
tant types of rare events are not accessible to any existing 
enhanced-sampling techniques, with the exception of the 
FFS-inspired method of Berryman and Schilling 1 ^ which 
relies on mapping the systems dynamics onto a time- 
inhomogeneous Markov process. The noise-sampling 
method of Crooks and Chandler^ allows sampling of 
transition paths in non-stationary systems but cannot be 
used to compute time-dependent transition rates, since 
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Figure 1. Barrier crossing in stationary state. Equilibration 
within the metastable region A occurs within ta- After a 
mean waiting time tab transitions from A to B occur, with 
a typical duration tc- 



that would require knowledge of the initial phase-space 
density, which is unavailable in general. 

In this article, we describe a new method, Non- 
Stationary Forward Flux Sampling (NS-FFS), that al- 
lows efficient simulation of rare events in time-dependent 
stochastic dynamical systems. NS-FFS constitutes a 
time-dependent generalization of FFS, and is concep- 
tually straightforward and easy-to-implement. NS-FFS 
achieves uniform sampling of trajectories crossing a pre- 
defined region in time and phase spac e, by combining 
interfaces in phase space as used in pFg 12 l 17 J with a flat- 
histogram pruned-enriched Rosenbluth method originally 
developed for polymer simulation d 18 * 19 ! In NS-FFS, tra- 
jectories are branched (proliferated) or pruned (termi- 
nated) based on their progression towards the final state, 
using interfaces in phase space and time. The scheme 
can be employed to sample the time-dependent phase- 
space density and time-dependent crossing fluxes, with 
uniform relative error. It thereby gives access to time- 
dependent transition rate functions^, including their 
low-propensity tails. 

The article is structured as follows. In section |TT] we 
provide a theoretical background, contrasting the well- 
studied setting of stationary, Markovian barrier-crossing 
with more general time-dependent rare event problems. 
Section [TTTJpresents the NS-FFS algorithm, together with 
corresponding pseudo-code. The correctness and effi- 
ciency of the algorithm are demonstrated in Section IV 
using two simple examples: diffusive escape in a one- 
dimensional W-shaped potential, and time-dependent 
switching in a genetic toggle switch. We conclude by 
discussing the main features of the method, and possible 
extensions, in Section fVl 



II. TIME-DEPENDENT RARE EVENTS 

A. Transition rate constants for Markovian systems in 
stationary state 

We first discuss rare transitions occurring in a system 
in stationary state. They may be visualized in terms of a 



static free-energy landscape^, as shown in Fig. [T] A typ- 
ical trajectory starts inside the metastable state A and 
fully explores the basin within a time ta\ after a wait- 
ing time tab it makes a rapid transition (on a timescale 
tc) over a high free-energy barrier C into state B. The 
essential observation is that if the equilibration time ta 
is much shorter than the mean waiting time in the A 
basin tab , then in the regime ta , tc *C t <C tab , tran- 
sitions from A to B occur in a Markovian, memory less 
fashion, effectively starting from a stationary state within 
A. Since the system is still in A with probability ~ 1, 
switches also happen with a constant propensity, whose 
value equals the rate constant Jzab — t~^ (see alscSSl). 
Numerical techniques for simulating rare events in sta- 
tionary systems exploit this fact by generating biased 
ensembles of short transition paths of duration ta,tc, 
which nevertheless allow reliable estimates of the much 
longer waiting time tab 3> TA) t C- 

We briefly review how this works in FFSW Given a 
progress coordinate A which increases from A to B, one 
defines a set of interfaces at successive levels {A;}i<;<l- 
A sample of points at the initial interface Ai is generated 
by a quasi-stationary simulation within state A. These 
are then used to initialize a set of trajectories which are 
propagated to interface A2 , or terminated if they re-enter 
A, whichever happens first. This procedure is repeated, 
starting from A2 and stopping at A3 or A, and so on. By 
propagating trajectories in segments from one A-interface 
to the next, FFS capitalizes on the rare fluctuations to- 
wards the transition. The resulting transition paths are 
of length tc < t <C tab- This leads to efficiency gains 
which grow exponentially with the barrier height. 



B. Time-dependent rare transition events 

In this paper we are interested in situations where the 
metastable states A and B can still be identified, but 
transitions between them happen with time-dependent 
propensity. For example, let us suppose that the generic 
system illustrated in Fig. [T] and discussed above is ex- 
posed to weak external forcing with protocol <f)(t), < 
t < T . For a macroscopic description with a time res- 
olution coarser that ta, the system is macroscopically 
Markovian and one can still define a transition rate from 
A to B, but this now depends on time: k,AB — kAB(tf^. 

Transition events may then be 'uniformly rare' so 
that k~^g(t) T for all t and the survival probability 
Sa (t) — 1 up to time T. However, if the transition rate 
is high in a particular time window, then the survival 
probability ^(i) may drop significantly below 1 dur- 
ing the time interval (0, T) of interest (as in figs. 9|10 
discussed below), and has to be taken into account in 
computing kAB{t^^- One then needs to measure both 
the first-passage time distribution or flux qAB(t) from A 
to B and the survival probability SU(i), to extract the 
time-dependent rate fc^g(t) = qabM / 'S A(t)', see the ac- 
companying pape: 
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Alternatively, transitions from A to B may be time- 
dependent even in systems without external driving, due 
to "macroscopic memory", in which the system's dynam- 
ics evolves on a time scale t s i ow such that ta < t s \ ow < 
tab- In this case, relaxation within the A basin will 
no longer be effectively instantaneous, and for t < t s \ ow 
the exit propensity from A will depend on the history of 
the trajectory. While transitions from A to B cannot be 
described by a rate constant, one may be able to char- 
acterize such systems in terms of a rate kernel kAB{t\t'), 
which quantifies the propensity to switch from A to B for 
the first time at time t, given that the previous switch 
happened at t' < t (see^l for a detailed discussion) . To 
extract the rate kernel from a simulation, one needs to 
measure the probability to stay in A without interruption 
from t' to i; this requires a simulation over a time interval 
(0,T) where T > t s i ow . In some systems memory effects 
may be combined with external driving - for example 
barrier escape in an underdamped, driven system. 

In all of these scenarios, the system dynamics is non- 
trivial and interesting over a time window (0,T) (deter- 
mined by the external driving or internal memory) which 
is longer than the typical transition time tab- This fact 
makes it impossible to speed up the simulation by gen- 
erating only an ensemble of short trajectories of length 
> ta , Tc, as is done in the rare event techniques for sta- 
tionary systems discussed above. To capture the physical 
behavior of non-stationary systems, trajectories must ex- 
tend over the entire time window (0, T) of interest. 

Nevertheless, enhanced sampling is both useful and 
feasible for non-stationary systems. Clearly, if in the time 
window of interest < t < T an event occurs with low 
probability, then brute-force simulations will fail to gen- 
erate more than a few, if any, events on this time scale. 
The goal of an enhanced sampling method is therefore 
to generate an ensemble of trajectories of full length T 
which is biased towards the transition. By reweighting 
this ensemble one can compute properties such as the 
time-dependent probability density, the transition flux 
Tab (t) , the transition rate function Uab (t) or the transi- 
tion rate kernel kAB(t\t') for the original system, over the 
time window < t < T of interest. NS-FFS, presented 
below, is precisely such a technique: by proliferating tra- 
jectories that evolve towards the final state and termi- 
nating those that do not, it generates transition events 
in the time window of interest; moreover, the branch- 
ing/pruning strategy is such that the relative sampling 
error is uniform over the space-time region of interest. 



III. NON-STATIONARY FORWARD FLUX SAMPLING 

The aim of the NS-FFS method is to generate a biased 
set of trajectories which sample transitions from state A 
to state B, defined by a progress coordinate A, as a func- 
tion of time, for non-stationary stochastic systems. To 
bias the set of trajectories towards the transition, one 
would like the flux of trajectories in the biased ensemble 



to be uniform in A (as in methods like FFS); to sample ac- 
curately the time-dependent behavior of the system (i.e. 
to obtain good sampling of early as well as late transition 
events), one would also like the flux of trajectories to be 
uniform in time. NS-FFS achieves both of these objec- 
tives, by generating a set of trajectories that uniformly 
cover a specified region of interest TZ in time and in the 
progress coordinate. Rare excursions into low-probability 
regions within TZ are sampled with the same accuracy as 
common events, and early excursions are sampled with 
the same accuracy as late ones. 

The method is conceptually simple. First, one defines 
the region of A-time space TZ that is of interest. One then 
partitions the region TZ using a series of interfaces; the in- 
terfaces may be defined as a level set cither in the progress 
coordinate or in the time coordinate. Each interface is 
then partitioned into a set of bins; the bins are defined 
as intervals in the respective other coordinate (time for 
A interfaces and A for time interfaces). 

The simulations start by generating an ensemble of ini- 
tial conditions, for example by performing a brute-force 
simulation in the initial state A. One then proceeds by 
firing off a trajectory from a randomly chosen initial con- 
dition, and propagating it according to the given dynam- 
ical equations of the system. The trajectory is assigned a 
statistical weight, which is initially unity. Upon crossing 
one of the pre-defined interfaces, the trajectory may be 
branched (split into several 'child' trajectories, with new 
statistical weights), or pruned (terminated). Repeating 
this procedure recursively for all child trajectories, one 
generates a 'trajectory tree' which extends form to T. 
The whole procedure is then repeated by firing off a new 
trajectory from a randomly chosen initial condition. 

The probability that a trajectory is branched or pruned 
when crossing a given interface bin depends on a running 
histogram which monitors the weighted flux of trajecto- 
ries crossing that bin. The branching/pruning rule is set 
up such that trajectories which arrive at a bin which has 
previously been under-sampled arc likely to branch; those 
which arrive at a previously over-sampled bin are likely 
to be pruned. The algorithm successively approaches a 
steady state in which the numbers of trajectories passing 
through all bins on all interfaces are equal - i.e. uniform 
sampling is achieved both in time and in the progress 
coordinate. 

In this section, we discuss each of the three 
ingredients — stochastic branching and pruning; inter- 
faces in phase-time space; the flat-histogram rule — in 
more detail, relegating mathematical details to Appen- 
dices. We then present pseudo-code of the resulting NS- 
FFS algorithm. Finally, we comment on the main fea- 
tures of the method. Again, implementation details are 
given in the Appendix. 
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Figure 2. Stochastic branching/pruning. A branching/ 
pruning move with average child number n proliferates or 
terminates branches. If n 1, surviving branches are de- 
creased or increased in weight, respectively (in the diagram, 
line thickness represents statistical weight). The target child 
number n may be an arbitrary function of the chosen coor- 
dinates q, q , as long as the reweighting factor r is satisfies 
Eq.[TJ 

A. Stochastic branching and pruning 

A key element of the NS-FFS algorithm is the branch- 
ing and pruning of trajectories, which allows control of 
the trajectory density. In a branching move, indepen- 
dent copies of the trajectory are created with a common 
history up to time t, while in a pruning move, the tra- 
jectory is terminated. The essential observation here is 
that one is free to branch or prune trajectories at will, 
as long as the statistical weights of the child trajectories 
are adjusted (reweighted) appropriately. 

Suppose that at time t, a trajectory with statistical 
weight w is randomly branched into n — 1 . . . n max chil- 
dren with probability b(n), or pruned with probability 
b(n = 0) (Fig. [2j. Each child branch is assigned a new 
weight w' — r(n) x w. Clearly this branching/pruning 
move will be statistically unbiased only if the weight fac- 
tor r(n) is chosen correctly. 

A necessary and sufficient condition for the combina- 
tion of b and r to be correct is that weight be conserved 
on average over branching/ pruning outcomes. That 
is, we may choose the branch number distribution and 
reweighting factor at will as long as they satisfy 

(nr) — b(n)nr(n) = 1, where b(n) — 1. (1) 

n=0 n=0 

This holds under very general conditions, including non- 
stationary system dynamics with memory and depen- 
dence of b on arbitrary parameters, as shown in app. [X] 
Thus we are free to adjust b to yield a desired mean 
branch number (n) £ (0,n max ), and thereby enrich or 
dilute the density of sample paths based on any chosen 
criterion, as long as we also adapt r to satisfy Eq. [TJ 

B. Interfaces in phase-time space 

In order to extend the FFS algorithm to non-stationary 
systems, we trigger branching/pruning moves whenever 



a trajectory crosses an interface. Since NS-FFS deals 
with systems whose dynamics are intrinsically time- 
dependent, the region of interest 1Z = [Xi,Xl] x [*i>*j] 
is two dimensional. To achieve uniform sampling of tra- 
jectories in both A and time, the branching/pruning rule 
needs to depend on both these coordinates. To accom- 
modate this, we define a set of interfaces as level sets in 
one of the coordinates (A or t), and partition each inter- 
face into a set of bins, which are intervals along the other 
coordinate. The region of interest 1Z is thus covered by a 
two-dimensional grid of subdivisions (see Fig. [3| . The in- 
terfaces are used to trigger the branching/pruning moves; 
the bins are used to determine the target child number 
n for these moves. 

The most direct generalization of FFS arises when in- 
terfaces are placed at a set of A-levels, and subdivided 
into time-bins (Fig. [3^,). The interface-bin grid is then 
given by the sets Bu = {(x,t)\X(x,t) = A;,tj < t < t i+ i} 
where bin Bu refers to the i-th time bin on the l-th X- 
interface. This interface arrangement will be referred to 
as 'A-if (A-interfaces). One can then measure the to- 
tal probability weight that has passed through bin Bu, 
Hu = ^2a=i Wa wnere -Wit is t ne running number of tra- 
jectories with weights {uu a } that have reached Bu. The 
crossing flux per tree is 

ju = H H /S, (2) 

where S is the running number of sampled trees. In com- 
puting Hu, one can choose either to count crossings only 
in the forward direction (increasing A), or in both di- 
rections, depending on the system property of interest. 
When counting only forward crossings, the quantity (ju) 
measures the forward probability flux across Bu', formally 

(ju) = J^ +1 (5(X(t) - x^xeix^dt, (3) 

where 9 is the Heaviside step function, and X9(X) is the 
forward probability flux. Although ju is strictly speaking 
a (unitless) crossing probability, it is proportional to the 
forward probability flux averaged over the bin Bu, which 
justifies the name 'crossing flux'. Importantly, when an 
absorbing boundary condition is imposed at the last in- 
terface Xl which is located beyond the top of the barrier, 
then 

(ju) - j 1+1 qA B (t)dt (4) 

where qAB is the first-passage time probability density, 
or exit flux, into B. This makes the A-if setup naturally 
suited for estimating exit fluxes. 

Alternatively, one may interchange the roles of time 
and the progress coordinate, by placing interfaces at a set 
of time points U, and subdividing them into bins along 
the A— direction (Fig.^). We call this interface arrange- 
ment 't-if (time-interfaces). The interface-bin grid is now 
described by Bu = {(x,t)\Xi < X(x,t) < Xi + i,t = ti}; 




Figure 3. In the A-if setup (a), branching/pruning events are triggered by the crossing of levels {A;} of the progress coordinate. 
Since in this case, the interface crossing flux ju corresponds to the flux of probability along A, the A-if setup is naturally suited 
for measuring fluxes and time-dependent transition propensities qA\(t). In the t-if setup (b), branching/pruning is triggered 
at a set of times {ti}. Here, the interface crossing flux ju provides a measure of the local probability density, making the t-if 
setup naturally suited for sampling p(A,t). 



bin Bu refers to the Z-th A bin on the i-th time interface. 
Crossing fluxes are still defined according to Eq. [2] but 
these have a different meaning in the t-if setup: ju sim- 
ply estimates the probability to find the system in the 
interval (A;, Aj+i) at time t^. 



{hi) = 



(5(X(U) - A))dA = 



p(A,ti)dA. (5) 



Thus the t-if setup lends itself naturally to estimating 
densities or potentials of mean force. 

In both the A-if and t-if setups, the branching/pruning 
rules are set up to ensure uniform sampling of ju- The 
A-if setup leads to uniform sampling of (forward) fluxes, 
while the t-if setup implements uniform sampling of 
phase space densities. We stress however that either 
setup could be used to measure either quantity. The 
relative efficiencies of the two methods depend on the 
quantity used for biasing but also on more technical as- 
pects, as discussed below. 



C. Sampling with uniform error 



trajectories that cross the bin, with an extra contribu- 
tion arising from the spread in their weights (Sw 2 ) /(w a ) 2 . 
Thus, to obtain a uniform relative error in ju (for a given 
total computational cost oc ^A^j), the branching rule 
needs to equalize the number of trajectories reaching each 
bin, while keeping the distribution of trajectory weights 
within each bin sharply peaked. In NS-FFS, these re- 
quirements are met by using a somewhat simplified ver- 
sion of the flatPERM rule™ 

Algorithm 1 Branching rule for bin Bu 

1. Calculate the target child number as n = w a /ju where 
w a is the weight of the incoming trajectory, and ju is 
the current flux estimate. 

2. Set the child number probabilities 

b(n) = 5 n [fi](n - [h\) + S n LR j ( \n] - n), 

where [■] and |_-J denote the ceiling and floor functions, 
respectively. 

3. Draw a child number n £ { [n\ , \n\ } from b. If n > 0, 
set all child branch weights to w' «— ju. 



The final ingredient in the algorithm is the rule for 
setting the child number probability b(n) for branching/ 
pruning moves. In NS-FFS, a target child number is set 
for each bin Bu, depending on the statistics of previous 
crossings of that bin; the goal of the branching rule is to 
sample the crossing flux ju through each of the bins Bu 
with uniform relative error. 

The relative error in ju in an NS-FFS run may be 
approximated as (see app. IbI 



a N 



(Ml 

(Ju) 2 ~ (Nu) 



1 



a v 



(6) 



Here w a denotes the (stochastic) weight of a trajectory 
reaching Bu, and averages and variances refer to an en- 
semble of NS-FFS runs with S trees. The constants 
aN,ce w approach unity in the 'ideal' case where trajec- 
tories which reach Bu are uncorrelated. This expression 
shows that the error is controlled by the number Nu of 



It is easily verified that the number of children is on 
average (n) = n = w a /ju- If n < 1, pruning occurs with 
probability l — h. The weight w' of the children (if any) is 
given by the parent weight w a multiplied by the reweight- 
ing factor r — w'/w a — l/n, which is independent of n. 
It is easy to see that the condition for unbiased statistics, 
Eq. [I] is satisfied by this branching/pruning rule. 

The branching rule, algorithm [T] produces a uniform 
error in the crossing flux estimate ju, since it tends to 
equalize counts between bins and minimize the weight 
variance within a bin. To see this, consider an NS-FFS 
simulation which is in steady state, i.e. after the crossing 
flux estimates ju have converged to their average values 
(ju)- In this situation the branching rule assigns each 
trajectory leaving bin Bu the 'perfect' weight w — > = 
(ju). Since all child trajectories are assigned this same 
weight, indeed the variance of trajectory weights leaving 
Bu is minimized (in the ideal case it is zero). We also 
note that the weight wff is equal to the system's intrinsic 
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probability to cross bin Bu. It follows that on average 
exactly one trajectory per tree (with weight wff) will 
emanate from bin Bu - or equivalently Nu/S x ju — > 1 X 
(ju) as the simulation converges. Thus NS-FFS achieves 
uniform sampling of the region TZ, with on average an 
equal number of trajectories crossing each bin, on each 
interface. 



D. The NS-FFS algorithm 

Combining the above ingredients, we now describe the 
full NS-FFS algorithm. To set up the simulation, one 
needs to generate initial configurations at time t — 0, ac- 
cording to an initial distribution p(xq) (which need not 
be known explicitly as a function of xq)', these could sim- 
ply be generated via a brute-force simulation in the A 
state. Next, one identifies a suitable progress coordinate 
A, and places a set of interfaces, subdivided into bins 
{Bu}, over the region TZ of phase-time space of inter- 
est, according to either the A-if or the i-if setup. One 
also specifies a dynamic range for re weighting, by setting 
minimum and maximum trajectory weights (w m i n , u> ma x) 
(these enhance convergence, see below). The simulation 
then proceeds according to algorithm [2j 

Algorithm 2 NS-FFS 

Init. Define a histogram of total crossing weights Hu and set 
Hu -v- for all I, i. Set a tree counter S <— 0. Define a 
queue of pending trajectories G and set G <— {}. 

Run Iterate the following steps, until the desired accuracy is 
reached: 

1. Start a new trajectory at t = 0, from a new initial 
state xq. Insert the trajectory into G, and as- 
sign an initial weight w <— 1. Increment the tree 
counter, S <— S + 1. 

2. While G is non-empty, iterate: 

(a) Pick and remove a trajectory from G. 

(b) Propagate the trajectory forward in time 
while recording any observables of interest, 
weighted with w, until either the final time 
T is reached, or an interface is crossed. In 
the latter case: 

i. Determine the bin Bu which was crossed 
and increment Hu by the current trajec- 
tory weight w. 

ii. If w £ (liimin, iBmai), carry out a branch- 
ing/pruning move: draw a child number 
n and set the child weight w' according 
to algorithm [T] Otherwise, set n <s— 1 
and w <s— w. 

iii. Generate n child trajectories and insert 
them into the queue G for further prop- 
agation. 



When setting up the region of interest TZ, clearly the 
transition region should be covered to enhance transition 



paths. It is equally important to let TZ extend well into 
the metastable basins; this allows for pruning of trajec- 
tories which would otherwise accumulate in the basins, 
degrading performance. Note also that by default, trajec- 
tories are not terminated when they leave TZ before they 
reach the final time T, ensuring that they may re-enter TZ 
and contribute at a later time. It is of course possible to 
explicitly add absorbing boundaries, which may increase 
performance, in cases where later reentry is not required. 

The weight limits (u> m i n , w max ) which appear in Alg.[2] 
are not strictly necessary for the existence of a steady 
state with uniform sampling of the crossing flux ju. They 
do, however, greatly enhance convergence towards it. In 
the initial phase of an NS-FFS simulation, the weight 
histogram Hu is sparsely populated such that the flux 
estimates ju are subject to large fluctuations. This can 
result in avalanches of correlated low- weight paths which 
degrade performance. For the method to be useful, it 
is necessary to have an effective way of controlling these 
bursts. Among a number of possible remedies including 
negative feedback control of tree size, branching thresh- 
olds, or explicit flat-histogram branching based on num- 
ber densities^, we found weight limits to be particu- 
larly simple and very effective. In practice, the reason- 
able rule-of-thumb to choose w m i n < xmni t i{{ju)} and 
w) max > 1 was found to work well. 

The output of algorithm [2] consists of weighted trees of 
trajectories in which all trajectory segments end either 
at an interface (branching/pruning point) or at the fi- 
nal time T (completion) (Note that a simulation may 
be stopped only after a full tree is finished). Trees 
may be generated depth-first (children before sisters), 
or breadth-first (sisters before children), depending on 
whether the queue G in Alg. [2] is of the last-in-first-out 
or first-in-first-out type. We found no significant differ- 
ence in performance between depth-first and breadth- 
first traversal, in contrast to other reports for the case 
of PER]Vp2l. Trajectory trees may also be written to disk 
in a recursive data structure for offline analysis. 



IV. APPLICATIONS 

A. Crossing of a linear barrier 

As a simple model for rare barrier crossing events, 
we consider overdamped Brownian motion in a linear 
double-ramp potential U = — a\x\ with boundaries at 
x = ±1 and barrier height a > 0, see Fig. [4] We consider 
two systems, which differ in their boundary conditions: 
one with reflecting boundaries at both x = — 1 and x = 1 , 
meaning that probability is conserved; another one with 
a reflecting boundary condition at x — — 1 and an absorb- 
ing one at x = 1. The two systems will be referred to as 
the reflecting/reflecting or reflecting/absorbing systems, 
respectively. This choice of model potential allows for 
an exact calculation of the Green's function, so that the 
correctness of our simulations can be assessed directly. 
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Figure 4. Linear ramp potential U (left) and the first 500 
branches of path trees generated during NS-FFS simulations 
using the t-if setup (top, reflecting boundaries) and the A— if 
setup (bottom, reflecting/absorbing boundaries). Although 
the total simulated time that is shown, is T sim < 6 <C T rxn , 
crossing events are successfully generated. Bins Bu are de- 
picted in blue, branch weight is indicated by line shading and 
branch points are shown as circles. 



Particles are injected atx=— 1,£=0 and diffuse 
according to the overdamped Langevin equation 



-DdJJ + Z. 



(7) 



We use thermal units where ksT = 1, so that £ repre- 
sents Gaussian white noise with (£(£)) = 0, (£(£)£(£')) = 
2D6(t — If). In these units, the diffusion constant D is 
equal to the mobility coefficient, and a has units of in- 
verse length. 

In this barrier escape problem both the time scale ta 
for equilibration in region A for x > — 1 and the crossing 
time tc, are much faster than the waiting time between 
crossings tab = t ba = 2t ix11 where T rxn is the global 
relaxation time of the system. The Laplace-transformed 
probability density p(x, s) = / °° e~ st p(x, t)dt can be cal- 
culated exactly using standard methods, see app.|D] lead- 
ing to explicit expressions for the time scales 



ta 



a 2 D 



aD 



and T r , 



2e a 



(8) 



In the present examples, the barrier height and diffu- 
sion constant were set to a — 15 and D = 1, respec- 
tively, such that (t a , t c , r rx „) = (0.02, 0.13, 2.9 x 10 4 ), re- 
spectively. Moreover, the full probability density p(x,t), 
and for the reflecting/absorbing case, the exit propensity 
j(x,t) 1 could be obtained as functions of time using an 
efficient numerical contour integration method^. 

Fig. [4] depicts the potential, the interface bins, and 
partial trajectory trees taken from NS-FFS simulations 



of this system from t = to t = T= 1. The region of 
interest was taken to be 1Z — [—1,1] x [0, 1] in (a;, f); the 
progress coordinate A was defined trivially as A = x. 

We first study the probability density p(X, t) in the 
reflecting/reflecting system, which relaxes towards the 



Boltzmann distribution 



over times longer than 



r rxn . The system was simulated using the t-ii setup which 
is most natural for measuring p(X,t), since the branch- 
ing factors are controlled by the local density. / = 199 
interfaces were placed across 7Z, at regular time intervals 
(from t = to 1), and partitioned into L = 40 equal-sized 
bins each (from A = —1 to 1). In Fig. [5 we compare 
p(X, t) as obtained as via exact analytical calculation, 
brute-force and NS-FFS simulation. While the proba- 
bility density p varies over 8 decades within the region 
1Z, the sampling density in NS-FFS is constant within 
a factor of 5 (Fig. |5j;). After reweighting, p is correctly 
reproduced throughout 1Z (even where it is very small; 
Fig. [5]i). This was achieved within a simulation time 
that would generate only a handful of transition paths in 
a brute force simulation (Fig. [EJj) . The region around the 
cusp of the potential is somewhat under-sampled. This 
is due to the fact that the force in our model is discontin- 
uous at x — 0; to achieve complete sampling uniformity 
in this region, one would require a bin width on the order 
of the length scale of variation of the potential. 

We next consider the time-dependent exit flux through 
the absorbing boundary in the reflecting/absorbing sys- 
tem. For this NS-FFS calculation, we used the A-if setup; 
this setup is natural since here the sampling bias is based 
on crossing fluxes over the A-interfaces, and the last cross- 
ing flux coincides with the observable of interest. L = 19 
A-interfaces were placed at regular intervals and parti- 
tioned into I = 50 time-bins each; only crossings in the 
positive positive A-direction counted towards Hu. Fig. [6] 
shows the probability density p(X, t) computed using the 
exact solution (a) and NS-FFS (c) as well as the un- 
weighted crossing fluxes and the reweighted exit flux (d) 
for this system. 

As expected, the number fluxes of trajectories ema- 
nating from each bin in the positive lambda direction 
(d, right axis) superimpose in a narrow band, confirming 
that NS-FFS indeed produces uniform sampling across 
1Z. Only the earliest bins at the farthest interfaces are 
visited less often, since the system dynamics does not 
allow them to be reached in the required time with suf- 
ficient probability. 

In contrast, the number density of trajectories over 
X,t (b) is roughly uniform over 7Z, but does exhibit a 
systematic bias favoring the A state. This can be un- 
derstood as follows. The branching rule in the chosen 
A-if setup biases towards uniform number flux in the for- 
ward A-direction, not uniform density. While in the re- 
gion x > trajectories spontaneously move in positive 
flux direction, in the region x < 0, trajectories tend to 
drift downhill, against the positive flux direction. In this 
region, the A-if NS-FFS simulation maintains a uniform 
population of uphill trajectories by proliferating. The to- 



8 




Figure 5. Time-dependent density p(A, t) for the linear ramp 
potential with initial condition p(A,0) = S(X + 1), with re- 
flecting boundary conditions, (a) numerically exact solution; 

(b) brute force sampling. In the t-if NS-FFS run, raw counts 

(c) cover 1Z almost uniformly while weighted counts (d) re- 
produce the exact density. Taking the B state boundary to 
be qs = 0.5, the occupation of the B state (d) fits a linear 
growth model with slope /cab = (3.4f ± .03) x 10~ 5 and delay 
0.099 ± .003 (blue, as extracted from NS-FFS; black, exact 
solution). Total simulated time was 10 ~ 3.4r rxn . Counts 
refer to the histogram bins of size (Ax, At) = (.05, .01) used 
in this figure. 



tal (uphill+downhill) trajectory density is thus increased 
in the region x < 0. 

This observation clearly demonstrates the difference 
in sampling biases between the t-if and A-if setups: the 
former generates a uniform density while the latter gen- 
erates uniform fluxes in the A-direction. Nevertheless, 
as Fig. [6]d shows, in practice the methods provide al- 
most uniform sampling of both forward flux and total 
number density, so that it is certainly possible to sam- 
ple fluxes using the i-if setup and densities with the A-if 
setup, without dramatic loss in performance. 




Figure 6. Time-dependent probability density (as in Fig. [5]) 
and exit flux, but for reflecting/absorbing boundary condi- 
tions. The exact density p is shown in (a). Counts generated 
in a A-if NS-FFS run show roughly uniform sampling (b, see 
text). Weighted counts (c) reproduce the exact density over 
the full dynamic range. Unweighted number fluxes of trajec- 
tories emanating from each bin (d, right axis) for all 19 A- 
interfaces (gray value increasing with A) collapse, indicating 
uniform number flux. As a consequence the exit probability 
flux j over x = 1 (d, left axis) is sampled uniformly includ- 
ing its low-probability onset (black, exact solution; blue, NS- 
FFS). An estimate of the stationary exit flux from these data 
for t > .24 gives j = 3.424(±.02) x 10~ 5 ; the exact value is 



Joo — t ab 



= 3.44 x 10" 



In this context it is worth noting that depending on 
the observable to be estimated, there is the additional 
freedom of setting up the simulation to measure cither 
density or exit flux. For instance, the slope in Fig. |5je) 
and the plateau value of the exit flux in Fig.|6jd) coincide, 
even though trajectories may re-cross the barrier from B 
to A in the reflecting/reflecting system, but not in the 
reflecting/absorbing system. This is in agreement with 
the exact expressions in the limit t <C r rxn ; indeed, as 
long as t <C r rxn , the occupation of the B state is negligi- 
bly small so that back-crossings occur with a probability 
even much smaller than forward crossings. Effectively, 
the B state initially appears as absorbing even in the 
purely reflecting system. Thus, either simulation may be 
used to measure the rate constant for this system. 

The relative errors in the estimated probability den- 
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Figure 7. Relative errors corresponding to the probability 
densities shown in figs. [5] and [6] for the t-if and A-if (top and 
bottom, respectively). Relative errors are computed as abso- 
lute differences between sampled and exact densities, normal- 
ized by the exact density. 



sity \Ap\/p= \p sim (X,t) -pexact(A,t)|/pexact('M) shown 
in Fig. [7] further illustrate the uniform sampling over 
TZ generated by NS-FFS. The relative error scatters uni- 
formly over TZ and in particular, does not scale with p -1 / 2 
as would be the case for brute-force sampling. Larger er- 
rors remain only in the fringes of the accessible region, 
where fewer trajectories are sampled. The residual stripe 
pattern in the t-if case carries the signature of correlated 
trajectories originating at the cusp of the potential; the 
undersampling right at the cusp which causes this can 
be considered a pathological feature of the force discon- 
tinuity in our model. Notice that although the A-if setup 
equalizes positive fluxes rather than number density, the 
probability density of trajectories nevertheless exhibits 
uniform error in this example, despite the weak asym- 
metry across TZ in the number of trajectories sampled, 
visible in Fig. [BJd. 



B. Genetic toggle switch 



As an intrinsically non-equilibrium example, we next 
consider a bistable gene regulatory network which can 
be seen as a simplified version of the A-phage genetic 
switch. This 'toggle switch' consists of two genes that 
mutually repress each other. In the 'exclusive' variant 
considered here, the two genes, which produce proteins 
A and B respectively, share a common DNA operator 
region O, such that when the dimer A2 is bound to the 
operator, protein B cannot be produce d, and vice versa. 
This model is discussed in detail in ref d 24 * 25 l 

In this simplified model, production of proteins is rep- 
resented by a Poisson process. The model consists of the 



symmetric reaction set: 

A' 



o 

A A 

fef 



A + A 



O + A 



A 2 
OA 2 



O 
B 

B + B 



A 
fet 



0+A 2 ^ OA 2 + B 2 ^ 

fe"off fc ff 

OA 2 4 OA 2 + A OB 2 A 



O + B 



B 2 
OB 2 
OB 2 + B. 



(9) 



The fact that only protein dimers may bind to the op- 
erator site, and the fact that the state OA 2 B 2 is dis- 
allowed, together make this system a robust bistable 
switch 24 . Each metastable state is characterized by an 
abundance of only one species, and transitions between 
these states occur on a much longer timescale than re- 
laxation within them. We use the same rate constants as 
irP^ \i = fc/4, fcf = fcb = 5fc, k on = 5k and k g = k, and 
we measure time in units of k^ 1 . 

A natural progress coordinate for this non-equilibrium 
system is given by the difference in total monomer num- 
bers, 

A = n B + 2n B 2 + 2n B2 - (n A + 2n A 2 + 2n A2 ). 

We are interested in a region TZ = {(A, i) e (—40,40) x 
(0, 10 3 )} which spans both metastable states and the 
transition region. Using the t-if setup, we define / = 500 
equidistant time interfaces. L — 16 A-bins were defined 
with boundaries at A = ±{40, 24, 22, 18, 15, 12, 9, 4} 
(These are the interface locations used irP^l, augmented 
by bins in the basins A,B.) We take the metastable 
states as A = {(x,t)\X < Xa} and B = {{x,t)\X > Xb} 
where — = Xb = 24. 



1. Unbiased relaxation 

Fig. [8] shows the result of an NS-FFS simulation of this 
model genetic switch from t = to T = 10 3 , with initial 
molecule numbers fixed to except n A = n^ 2 = 10, so 
that A(0) = —30 € A. As in the previous one-dimensional 
example, the region of interest is sampled approximately 
uniformly in NS-FFS. Measuring the occupancy of the B 
state (fiB(t)) = (9(X(t) — A_b)), and fitting to a delayed 
linear rise (t — ri ag )/rAB, we obtain a lag time ri ag = 
(129 ± 5), and recover a waiting time for barrier crossing 
t ab = (1.07±.01) x 10 6 in accordance with the previously 
measured value l/fc A B = (1-06 ± .02) x 10 6 ^ 



2. Response to time-dependent forcing 

We now consider the reaction of the toggle switch to a 
time-dependent external bias. This case is inspired by the 
phage-A switch in the bacterium Escherichia coli, where 
an increase in intracellular RecA concentration triggers 
the transition from the lysogenic to the lytic phase of 
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Figure 8. Bin counts (a), probability density p(\,t) (b), 
and cumulative crossing probability (/ib(£)) (c) for the toggle 
switch, simulated using NS-FFS with the i-if setup. A linear 
fit for t > 250 is shown in orange. Error bars were generated 
by bootstrap resampling from 25 independent simulation runs 
with total simulated time 10 6 each. 



the virus life cycled As a simplified model for the ac- 
tion of RecA we introduce a species R which degrades A 
monomers: 



R^ 



A 



R ^> R 



(10) 



The degradation of A by R forces the switch towards the 
B state; thus R can be regarded as an 'external force' 
acting on the switch, whose strength can be measured 
by the steady-state bias 7n£P, where n£? = /cr/Vr is the 
number of molecules of R in steady state. Relaxation 
of 7JR towards is exponential with a relaxation time 

Mr 1 - 

First, we initialize the switch in state A with ur = 
copies of R and switch on the production of R. The 
steady state bias is chosen such that in steady state 
the switch is fully driven to the B state. The switch 
then flips from A to B with a distribution of switching 
times. Switching events result from favorable fluctua- 
tions in the copy numbers of the molecules that consti- 
tute the switch (eqs. |9j. These are partly due to the 
intrinsic stochastic nature of the switch, and partly in- 
duced by (extrinsic) fluctuations in the number of bias- 
ing molecules R (eqs. 
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The distribution of switch- 
ing times thus reflects both intrinsic noise of the switch 
and extrinsic noise originating from fluctuations in the 
number of R molecules. We investigated these effects 
by using NS-FFS to obtain the switching dynamics as a 
function of the level of noise in the bias. To modulate 



the latter, we varied the equilibrium copy number of R 
between — 1 and = 100, while keeping the av- 
erage bias 7n£? and bias time constant /i^ 1 = 500 fixed. 
Therefore, on a mean-field level, all bias protocols were 
kept the same. However, the individual bias trajecto- 
ries n-R(t) are markedly different: At — 100, each 
trajectory nn(t) exhibits a nearly deterministic and ex- 
ponential rise in time, while at — 1, each individual 
trajectory 71r(£) is a single off-on event, which is expo- 
nentially distributed in time. Fig. [9^, shows the mean 
exponential rise of the bias protocol and the level of fluc- 
tuations around it. Note that in a linear system where 
the state occupancies are linear functions of the external 
forcing history, this variation of parameters would lead 
to a probability of having switched after the pulse which 
is independent of the pulse duration /j,^ 1 . 

Fig. |9jb-e) illustrates the switch response to a bias at 
various noise levels. For a nearly deterministic bias, with 
ng 3 = 100, the switch response is characterized by a grad- 
ual increase of A towards a threshold, followed by a tran- 
sition over the threshold (Fig. [9Jd) . The transition times 
have a relatively narrow distribution around t = 500. 
Since further increasing the expression level does not 
sharpen the switching time distribution (not shown), the 
result shown in panel b corresponds to the intrinsic limit 
in precision of the toggle switch at the given rate of bi- 
asing /!R. 

As the driving becomes noisier (Fig. [9j:,d), the switch 
response exhibits a broader distribution of switching 
times, with both early and late crossings, and in addition, 
a weak new metastable state develops at the transition 
state (d) . Fig. [9^ summarizes these results, showing the 
probability that the switch has flipped, as a function of 



time, for several different values of 



It is clear that 



noise in the driving force has an important effect on the 
switching trajectories. 

The qualitative change in the switch flipping trajecto- 
ries for low values of n^P can be understood by a picture 
in which the switch dynamics are enslaved to fluctuations 
in R. In the extreme case n^P = 1 (Fig. |9jl), individual 
bias trajectories switch to full bias strength suddenly, at 
random times. The switch then responds to the strong 
bias by rapid nipping to the B state in a stereotyped 
way; the time course of the switching event (not shown) 
includes a pause at the threshold which is responsible 
for the zone of higher density around A = (This pause 
originates from the fact that for n£? = 1, the A molecules 
are very rapidly degraded when the single R molecule be- 
comes present, while it still takes time to produce the B 
molecules, which ultimately flip the switch). 

We now consider the switch response to transient 
pulses of biasing molecules R. In these simulations, af- 
ter the system reaches a quasi-stationary state in re- 
gion A, ?ip Ulsc = 10 2 biasing molecules are flushed in 
instantaneously. We set A;r = 0, so ur then decays 
to over a pulse duration /i^ 1 . During the pulse, 
the switch is biased towards the B state. To isolate 
the effect of different pulse durations on the switch re- 
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Figure 9. Response of the toggle switch to a bias of R 
molecules, whose production begins at time t = (Eq. [10 1. 



The mean bias strength approaches 7?!^ = k (a, lines) but 
with different amounts of noise, corresponding to different 
choices for n£? (a, shaded areas lie between the 5th and 95th 
percentiles). The switch response becomes more random as 
the noise in the bias increases (n^ decreases; b-d, respec- 
tively), and the transition time distribution widens (e, error 
bars indicate the 5th and 95th percentiles for the B state 
occupation (/is(t))). 



sponse, we adjust 7 such that the integrated bias strength 
(/ 7 n R(£)di) = 1 remains fixed, while changing the value 



of 



/'R 



Or 



100, 1,0.1). This leads to 7 = MrA% 



pulse 



Fig. [10] shows that the toggle switch possesses an op- 
timal pulse duration, even though the integrated bias 
strength remains constant. For moderately long pulses, 
the occupation of the B-state after t = 500, increases 
with decreasing pulse duration. This can be understood 
in terms of a non-linear threshold behavior of the switch. 
As the pulse duration fj,^ 1 is decreased, the initial bias 

strength 7^r u1sc = /Jr increases, enhancing the switching 
probability. However, as the pulses become even shorter 
the switching probability decreases again. This decrease 
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Figure 10. NS-FFS results for the toggle switch, under the 
influence of pulses of biasing molecules R of varying duration 
/i^ but constant total efficiency {J ynn(t)dt) = 1. The time- 
dependent bias {^yn-R.) is shown in (a). The response of the 
system to pulses of durations fi^ = {100, 1,0.1} is shown in 
(b,c,d), respectively. The crossing probability, taken to be the 
B state occupation {hs{t = 500)} exhibits a maximum at a 
pulse duration around fi^ = 1 (e) 



is a dynamical effect: since the switch cannot respond to 
changes in riR which occur on a timescale shorter than 
its own intrinsic kinetic time scale it acts as a low- 
pass filter. In this sense the toggle switch can be said 
to be robust against both strong transient perturbations 
and persistent weak perturbations. These results clearly 
show that the response of a genetic switch to a pertur- 
bation (or signal) is a dynamic property which depends 
not only on the (integrated or peak) pulse strength of the 
perturbation but also on its pulse shape. 



V. DISCUSSION 

A. General features 

The NS-FFS scheme has a number of characteristic 
features. First, like FFS, NS-FFS does not perturb the 
given dynamical equations of the system (i.e. no bias- 
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ing force is applied). Instead, NS-FFS generates a bi- 
ased ensemble of unbiased trajectories by proliferating 
those that move in a preferred direction and terminating 
those that do not, with appropriate reweighting. This 
means that NS-FFS, like FFS, is suitable for systems for 
which the dynamical equations cannot easily be biased 
and reweighted (e.g. because they do not obey detailed 
balance). It also makes NS-FFS highly suitable for im- 
plementation as a wrapper around existing simulation 
code. 

Second, NS-FFS generates trajectories with a complete 
history from t = on. This means that no assump- 
tion of memory loss is made when a trajectory passes 
between interf aces , as in some other rare event simu- 
lation methods^'. Perhaps more significantly, NS-FFS 
also does not assume loss of memory on reentry to the 
A state. Most existing rare-event m ethod s for stationary 
systems, such as TPS, TIS, and FFS^^, rely on the as- 
sumption that trajectories which re-enter A equilibrate 
rapidly and can be treated as new, independent trajecto- 
ries when they eventually re-exit A. For most stationary 
systems where ta <C tab this is a reasonable assumption, 
but for non-stationary systems, such as those with non- 
Markovian or time-inhomogeneous macroscopic switch- 
ing dynamics^, correlated entrance and re-exit from a 
basin can make a significant contribution to the time- 
dependent quantity of interest. Thus it is an essential 
feature of NS-FFS that memory loss is not assumed for 
t > 0, even if the system re-enters the A basing. 

Third, once the NS-FFS algorithm has reached its 
steady state, each interface bin emits one trajectory per 
started tree on average, so that trajectories are sampled 
uniformly over the range TZ (see figs.[5][6j [8]). Some fluc- 
tuations around this uniform sampling average arc toler- 
ated in exchange for narrow weight distributions at each 
bin. These features are a direct consequence of the basic 
flatPERM branching rule^, which does not implement 
a negative feedback on (unweighted) trajectory numbers. 
NS-FFS is thus a 'weak' flat-histogram method. 

Fourth, we note that the effectiveness of a simulation 
scheme depends not only on its ability to generate many 
samples but also on the independence of these samples. 
Clearly, in NS-FFS, after a branching event, child trajec- 
tories remain correlated for a certain time. This suggests 
that one should allow further branching of the children 
only after a refractory time of the order of the typical 
decorrelation time. We did not observe this modification 
to produce any significant improvement for the systems 
studied here. This is presumably because these systems 
are sufficiently stochastic that branched trajectories any- 
way decorrelate rapidly between interfaces. In contrast, 
it turned out to be crucial to control the exponential 
growth of trajectory trees in the initial phase of a simula- 
tion since these generate highly correlated samples, which 
delay convergence of the crossing weight histogram. This 
is simply and efficiently accomplished by using weight 
limits, as described above. 



B. Progress coordinate vs. time based branching 



The A-if and f-if interface setups (sec. Ill F3 1 differ in 
that the former equalizes number fluxes of trajectories 
in A-direction across the region TZ while the latter equal- 
izes their number density. Nevertheless, it is of course 
possible to use either scheme to measure any quantity 
in a given physical system; the choice of setup will af- 
fect only the efficiency of the calculation. We now briefly 
discuss how the observable of interest and the computa- 
tional overhead associated with branching can affect the 
choice of the most suitable setup. 

Target observable Suppose one wishes to measure the 
time-dependent propensity &ab(£) for exit over some fi- 
nal level Xl of the progress coordinate at the boundary 
of state B, over some time interval [ti,tj] (cf. Fig. [IJ. 
If re-crossings back from B can be safely neglected, we 
may place an absorbing boundary at Al- It is then nat- 
ural to use the A-if setup, placing A-interfaces at levels 
{Xi}i—i...L, and count forward crossings only. In the limit 
of slow escape over Xl such that the survival probability 
J p(X, t)d\ ~ 1 over the time of observation, the observ- 
able kAB coincides with the positive flux over the last 
interface A/j^l Since the relative error in the positive 
flux is equalized over all preceding interfaces, the A-if 
setup will generate uniform sampling for fcAs(i) at the 
final interface over the time interval of interest. 

Alternatively, one may be interested in a potential of 
mean force — logp(X,t) over a region TZ = [Xi,\l] x 
[ti , ti] in A— t space. In that case the t-i£ setup is the more 
natural choice, since it generates a uniform relative error 
in p(A, t) over all bins. One then obtains an estimate of 
— log p with uniform absolute error over the region TZ. 
In particular, saddle points and basins are sampled with 
equal frequency. 

Branching overhead The A-if and i-if setups differ 
in the relative overhead of coordinate evaluation and 
branching. Clearly, the detection of crossings requires 
the evaluation of A, which incurs some computational 
overhead o\ per evaluation. The branching move itself 
and the maintenance of the tree structure add a second 
type of overhead Of,, which we also take to be constant 
per branching event. 

In the i-if setup, A evaluation and branching are cou- 
pled and happen once every time-interface spacing tj; 
the overhead is (o\ + Ob)/T t per unit simulated time. In 
contrast, in the A-if setup, they are decoupled: interface 
crossings are checked at intervals t\ while the trajectory 
is being propagated, whereas trajectory branching hap- 
pens only if a crossing is detected. The A-evaluation in- 
terval t\ is an independent adjustable parameter (but 
must not be too long or the algorithm will fail to detect 
crossings). If interface crossings happen every n on av- 
erage, the overhead is 0\/t\ + o^jr^ per simulated time. 

If A evaluation is expensive (p\ ^> Ob), it is useful to 
increase t\ or r t , respectively, as far as possible without 
degrading the weight statistics. In the opposite case that 
A evaluation is cheap (o\ < oj), the A-if setup may be 
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advantageous, since this setup gives the option to check 
A often (t\ <C n), at small cost. Especially in systems 
where excursions towards the B state tend to be short- 
lived, it is advantageous to check for interface crossings 
significantly more often than they actually occur, since 
this increases the chance of detecting and capitalizing on 
short-lived forward excursions. The dependence of NS- 
FFS performance on parameters will be addressed more 
fully in a future publication. 



C. Variants and extensions of the algorithm 

Relation to stationary FFS If NS-FFS is used to sim- 
ulate a system which is in stationary state, early and 
late barrier crossings are equivalent; crossing statistics 
can then be improved by binning early and late crossings 
together. This can be achieved within the A-if setup: 
trajectories are started from a stationary distribution at 
t — 0, and a single time bin is defined (I = 1), lead- 
ing to quicker convergence of the, now one-dimensional, 
crossing weight histogram H — {-ff;i}i<i<L- This sta- 
tionary version of NS-FFS is similar but not identical 
to the 'branched growth' variant of conventional FFS 1 '. 
In branched-growth FFS, trees have a fixed number of 
children at each A interface crossing and are terminated 
exclusively at the next interface or when the system re- 
turns to A; the stationary NS-FFS can be seen as an 
adaptive generalization of this scheme. 

Multiple progress coordinates In problems with multi- 
ple alternative transition pathways, finding a single good 
progress coordinate can be challenging. Use of a poor 
progress coordinate fails to distinguish between trajecto- 
ries which are likely and unlikely to result in a transition, 
making successful biasing of transition paths impossible. 
To some degree, NS-FFS already alleviates this problem 
compared to TIS or FFS: the extra dimension of time 
acts as a second progress coordinate which allows us to 
discriminate between early and late crossings. In a class 
of systems which show distinct 'slow' and 'fast' pathways 
of a reaction between A and B, an NS-FFS simulation 
would be able to separately enhance these pathways. 

If different transition pathways share a common time 
scale then time as an additional progress coordinate is not 
useful in itself, but one may be able to find a small set of 
progress coordinates which successfully separate different 
pathways. In this case, one can make a straightforward 
generalization of the NS-FFS scheme to multiple progress 
coordinates. As in any multi-dimensional scheme, this 
will come at the cost of more bookkeeping and possibly 
slower convergence of the weight histogram H due to the 
increased number of bins, see app. [C] 

Parallel version The NS-FFS algorithm lends itself to 
a parallel implementation. Each trajectory can be simu- 
lated independently until an interface is crossed. At this 
point, the shared histogram Hu is read and updated. Af- 
ter branching, n parallel simulations for the children are 
spawned. The communication between simulation pro- 



cesses is restricted to histogram updates; depending on 
available bandwidth these updates could also be cached 
and applied in groups, without biasing the sampling. As 
the global histogram converges, updates become unneces- 
sary and the simulation gradually becomes trivially par- 
allel. 

Adaptive generalizations As shown in app. [A] branch- 
ing/pruning events may be introduced at will as long as 
weight is conserved on average (Eq. [TJ. This includes 
complete freedom of: adaptive updates of the bin bound- 
aries or the interface placement; inserting or removing 
interfaces; smoothing of bin counts within or across in- 
terfaces; or pre-filling of crossing histograms based on 
prior knowledge. 

All of these options should allow for further perfor- 
mance improvement in particular situations, to be ex- 
plored in future work. In particular, it is interesting to 
ask if an optimal interface and bin arrangement can be 
found iteratively, as has been proposed for FFS-^l A 
promising direction might be to monitor the local sam- 
pling noise and adapt the interface arrangement in re- 
sponse to it. 



VI. CONCLUSION 

In this article, we have introduced an enhanced sam- 
pling scheme, called NS-FFS, which is conceptually sim- 
ple allows the efficient sampling of rare events in non- 
Markovian and non-stationary systems. 

The NS-FFS algorithm builds on two widely used in- 
gredients: a flat-histogram branched growth algorithm 
closely related to PERMASED anc j the concept of phase- 
space interfaces^ to monitor progress towards a tran- 
sition. NS-FFS is a generalization of FFSP^, and is 
straightforward to implement, especially when one does 
not want to store the trajectories themselves. We have 
demonstrated the correctness of the method, and given 
several simple example applications which highlight both 
the effectiveness of the method and the relevance of in- 
trinsically time-dependent rare events. 

A host of physical, chemical and biophysical prob- 
lems are amenable to NS-FFS simulations. These in- 
clude the computation of time-dependent transition rates 
in systems with time-dependent external driving^ such 
as the signal- induced flipping of genetic switches studied 
here, crystal nucleation during a temperature quench or 
protein unfolding under force; non-exponential switch- 
ing time distributions in processes that can be coarse- 
grained as switches with me mory, such as the switching 
of the bacterial flagellar motoil^SEll and escape probabil- 
ities from a non-equilibrium distribution in a metastable 
initial state within a prescribed 'window of opportunity', 
like the flipping of genetic switches induced by transient 
pulses. 
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Appendix A: Reweighting criterion 

In this section we sketch a proof that the condition of 
weight conservation on average over branching outcomes, 
Eq. ^ is sufficient for unbiased sampling in a general 
branched-tree simulation scheme. This amounts essen- 
tially to careful bookkeeping. No assumptions are made 
about stationarity or loss of memory. 



1. Statement of the problem 

Consider a stochastic process x(t) which is started at 
t = with a value of x(0) drawn from some initial dis- 
tribution. The process is fully characterized by all of its 
m-point joint probability density functions (pdfs) 

p(x m ,t m ; 5Gl,tl) = (S(x tm - x m ) ■ ■ ■ 5(x tl - Xx)), 

(Al) 

for the trajectory to pass by the sequence of sample points 
at times < t\ < t 2 < ■ • ■ < t m . Here x tm 
is a shorthand for x(t m ). In a brute-force simulation, the 
joint pdf is estimated by an average over S independent 
runs, 

1 S 

p(x m ,t m ;...;xi,ti) ~ — }J(x s tm - x m ) ■ ■ ■ 5{x s ti -xi). 

s—l 

(A2) 

(we note that the right hand side of this equation is sin- 
gular; the approximate equality is implied when Eq. |A2| 
is integrated over finite regions). 

We now introduce a single branching event at an 
intermediate time t' , and let m! — max{m|i TO < 
t'y. Upon branching, n' statistically identical copies 
of the system with independent futures are gener- 
ated. That is, for a given history up to time 
if, each of the copies has the same conditional pdf 
P\x m , t m ; . . . ; x m >+i, t m i-^i\x m i , t m i ; . . . ; X\i t\) to visit 
future points of phase space, but future points of dif- 
ferent children are mutually independent. No Markov 
assumption about the system is being made. The task 
is now to assign correct weights to the child branches in 
order to guarantee unbiased sampling. 



2. Unbiased sampling with one branching event 

An obvious choice for the weights is to conserve the 
total in- and outgoing weight at the branching point and 
to treat children equally. This rule gives a weight 1 jvf 
for each child trajectory from if on. 

However, this rule disallows pruning, i.e. n'=0. To 
enable pruning, it is necessary to relax strict weight con- 
servation at the branch point. To do this, child numbers 
n' = 0,1,... n max are drawn at random with probabili- 
ties b(n'), X)n'=o^( n ') = 1- We then assign weights w' c 
to the child branches c = 1 . . . n' (if any). The expected 
total weight of trajectories passing through a sequence of 
points is then 

W(x m ,t m ; . . . ;xi,ti) 

n' 

= (XI W 'c S ( X tm - X rn)--- K X t m , +l ~ 2W + l) 
c=l 

xS(x tm , - aw) • • -S(x tl - xi)^ 

n' 

= (X W ( S ( x t m _ x m) ■ ■ ■ S{x fl - Xi)) 
c=l 

n 

= (^2 w 'cjPi x m,t m ; . . . ;xi,ii), (A3) 

c=l 

where mf is defined as above, sums running from n' — 1 
to nl = vanish by definition, and we used the fact that 
children are identical. Clearly the condition of weight 
conservation on average, 

Mi>0= 5> (n,) i>- (A4) 

c— 1 n' — l c—1 

is necessary and sufficient for unbiased sampling, since 
then W = p, i.e. the reweighting has corrected the bias. 

Since there is no reason to treat child branches differ- 
ently, we set w' c = w' = r(n)w where w = 1 is the parent 
weight. Eq. |A4| now reduces to 

^max 

1 = (nV(n')) = £ b(n')n'r(n'), (A5) 

n' = l 

Eq. |A5| constrains the choice weight factors r for a given 
branching distribution b. For instance, if n max = 2, and 
(6(n)) n= o,i,2 = (-5, .2, .3) then the choices {r{n)) n= x,2 = 
(0,5/3), (5,0) and (5/4,5/4), are all unbiased. 

Generalizing Eq. |A2| we can then estimate p from a 
branched simulation of S independent trees as 

piXjyi , t m , . . . , Xi , tij 
, S n' a 

8=1 C s — 1 

x 6{xl m , - x m >) • • • S(x a tl -xi) (A6) 
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The first line on the rhs contains sample points after 
branching into n' s children (if any), and the second line 
those before branching (if any) ; the weight factors r sat- 
isfy eq. |A5| Note that repeated simulations average over 
not only the system but also the branching randomness. 
For instance, the total weight of all trajectories at a given 
time is conserved only on average over branching out- 
comes. 

So far we have shown that a single stochastic branching 
move at time t' and re-weighted sampling according to 
Eq. |A6| does not introduce a bias if the child number 
probabilities b(n') and child weight factors r(n') obey the 
condition of weight conservation on average, Eq. |A5[ We 
note that to verify this numerically, one would need to 
generate many trees always with branching time t' and 
count joint hits of all bins around the points x\, . . . , x m 
at the times t%, . . . , t m ; if the final time is larger than 
the branching time, hits are re- weighted the factor r(n s ) 
appropriate for the respective branch number n s . 



Each new branching event with n l children adds a 
weight factor r(n l ), so that the instantaneous weight 
along a trajectory passing by the branching events 
(n', t'\ n 2 , i 2 ; . . . ; n k ,t k ) becomes 



w(t; n , t ; 



■n\t k )= [] » 

l:t l <t 



(A8) 



The system density estimate Eq. |A6| now takes the 
form of a hierarchical sum over all weighted branches 
existing at the final time. This is simplest to write 
down using a recursive definition. We index a particu- 
lar branch (trajectory segment) anywhere in the tree by 
the sequence of the child indices starting from the root, 
7 = (c, c', . . . , c k ), and let |7| = k + 1 denote its nesting 



depth. The zeroth child index c 



s is the tree 



index. The points along the trajectory from the root up 
to and including 7 are denoted xj . Let 



7T 7 (£r, 



;^i>*i) = n 8(3 



(A9) 



3. Adaptive branching probabilities 



if 7 has no children, and define recursively 



We now show that the branching probabilities may 
be adapted according to an arbitrary protocol. To see 
this, let b depend on an arbitrary parameter a such that 
^Zn'bain') = 1, and choose r a accordingly, such that 
~Yl, n ib a (n')n'r a (n') = 1 for all values of a. Then, fol- 
lowing the preceding discussion, a branching event is un- 
biased for any a. We may take a to be any random 
variable, depending on a set of bias control parameters /3 
via some density p{a\f}). Since Eq. A3 now reads 



W(x m ,t m ; . . ■;x 1 ,t 1 ) 

= J dap(a\/3) ^ b a (n')n'p(x m ,t m ; . . . ; xi,ti)r a (n') 



dap(a\/3)p(x m ,t m ; . . .;xi,ti) 

= p(x m ,t m ;...;x 1 ,t 1 ), (A7) 



unbiased sampling is still guaranteed. We thus have com- 
plete freedom to introduce dependency of the branching 
probabilities on arbitrary extra information, including 
the past or future system state, or the history of the 
simulation. 



4. Multiple branching events 

Finally, we generalize the arguments above to multiple 
independent and non-synchronized branching events on 
different branches of a tree. The idea is that introduc- 
ing a new branching event at a time t" on an existing 
child branch is bias-free, as long as the new b and r fulfill 
Eq. |A5| We then argue by induction over tree genera- 
tions. 



^ {Xm ) tm j • • • j 3Z\ , t\) ^ ^ 7T ^ ^ {x m , t m , . . . ,Xi,ti), 

' (A10) 
if 7 has n 7 children with weight factor r 7 . This recursion 
terminates since the simulation has a finite maximal nest- 
ing depth K. Note that 7r 7 (x m , t m ; . . . ; X\, tx) counts the 
weighted hits of 7 and all its descendants to the points 
(x m , t m ; . . . ; Xi, ti). The weight is effectively given by 
Eq.gg] 

All S simulated trees together can then be represented 
as the descendants of the empty path (). If we define 
Tq = l/S, tiq = S and |Q| = 0, the estimate for the 
joint m-point density is 

p(x mi t m ; . . .;xi,ti) ~ {Ti { \x m ,t m ; . . . ;xi,ti)). (All) 



Since Eq. |A9| coincides with the brute force estimate of 
the density, and Eq. |A10l is unbiased by the arguments in 
sec. |A 2| it follows by induction over K that the density 
estimate E q. |All|is u nbiased. 

Eqs. ( A9 AlO|Alf ) directly translate into a histogram- 
ming algorithm, which recursively parses a trajectory tree 
while binning weighted counts. Binning can be carried 
out on-line. This is straightforward if only the one-point 
density p(x, t) is desired; online updates of m-point den- 
sities would require m—1 nested inner loops over the tree 
for each simulation step. 

As a corollary, arbitrary m-point observables A can be 
estimated along the same lines. One just has to replace 
7r by n A in Eq. |A10| and|A"9|by 



\t m ) = A(x tm ,. . . ,x tl ), 



(A12) 



if 7 has no children. For instance, a two-point auto- 
correlation function would correspond to A(x t2 ,x tl ) = 
x t2 x tl - (x t2 )(x tl ). 
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Appendix B: Variance of weighted bin counts 

Here, we write the variance of the total weight W — 
^2a =1 w a accumulated in a given bin Bu in terms of the 
statistics of the trajectory weights w a and the number N 
of trajectories which has reached bin Bu. We recall the 
'law of total variance': 

(A 2 ) - (A) 2 = ((A 2 \C) - (A\C) 2 ) + (((A\C) 2 ) ((A\C)) 2 ) 



Here the inner expectation values are conditioned on 
some event C, and the outer expectations average over 
C. We will denote the conditional variance as (S X 2 \C) = 



Eq. 



Bl 



then be- 



(X 2 \C) - (X\C) 2 for any observable X. 
comes (5A 2 ) = ((SA 2 \C)) + (S(A\C) 2 ). 

Consider an idealized NS-FFS simulation in which tra- 
jectories arriving at Bu are uncorrected. Specifically, we 
assume that the incoming trajectory weights {w a } are 
mutually independent and independent of N, and that 
arrivals are a Poisson process. For the mean total cross- 
ing weight we then obtain (W) = (N)(w a ). The total 
weight variance is given by 



(6W 2 



((SW 2 \N)) + (S(W\N) 2 ) 
(N)(Sw 2 a ) + (w a ) 2 (SN 2 ) 
(N)((Sw 2 a ) + (w a ) 2 ); 



the relative variance of the collected weight becomes 



(SW 2 

(wy 



1 

W) 



i 



In the more general case, we now assume that correla- 
tions between branches increase the noise while respect- 
ing the same scaling with the count number, and write 



(SW 2 

(wy 



(N) 



1 



(B2) 



where a at > 1 if trajectories arrive in bunches and a w > 1 
if their weights are correlated. The crossing flux ju is 



estimated as Wu/S and thus has the same noise, Eq. B2 
as Wu. 

The noise in bin weights can be split up further. The 
incoming weights {w a } are distributed with a mean and 
variance which result from both the inter-bin variance be- 
tween starting bins £?/v and from the intra-bin variances 
of outgoing weights from within B^i* . We have 



(B3) 



(w a ) = ((w a \l'i'}), and using Eq.[BT 
(5w 2 a ) = ((6w 2 a \l'i')) + (6(w a \l'i') 2 y, 



here the first term is the mean intra-bin variance within 
originating bins, and the second term is the inter-bin vari- 
ance. Plugging in Eq. |B3| we can write the noise in W 
as 



(SW 2 ) 

~(wy 



W) 



(WTO 

(w a ) 2 



(S(w a \l'i'} 
(w a ) 2 



As the simulation progresses, crossing flux estimates 
converge, so that all trajectories leaving i?;v are even- 
tually assigned the same weight. The intra-bin term 
((8w 2 a \l'i')) / (w a ) 2 cx (Sjli,) /(ji>i') 2 thus vanishes as N -> 
oo. The inter-bin term (8(w a \l'i') 2 ) / (w a ) 2 reflects the 
non-uniform transition probabilities between bins and 
persists also in steady state. 

In order to balance noise contributions, it seems rea- 
sonable to choose bin size and interface spacing such that 



(Bl) a w {5{w a \l'i') ) — (w a ) m an equilibrated simulation 



Appendix C: Multi-dimensional progress coordinates 

A multi-coordinate NS-FFS simulation can be set 
up as follows. First, find K progress coordinates 
{^ k }k=i...K, with corresponding sets of levels {A§ < • • • < 
\ k Lk \k=\...K- Denote the interval (Af_ x , X k ) = Aj, . Define 

a subset of K' < K progress coordinates {X k ' }k'=i...K' ■ 
Only the first K' sets of interfaces will trigger branching 
events; for these interfaces, define bins: 



B 



k' 

tli2-..ifc' ...1} 



{(x,t)\X k (x,t) = Af , and \ k (x,t) G Af k for k + k'} 

(CI) 

Then, proceed as before: branching moves are triggered 
on the first K' interfaces, based on the corresponding 
crossing weights H^ h ly 1r . 

In this setting, time is treated as another progress coor- 
dinate. To recover the NS-FFS setups discussed in sec. |III| 
in this setting, let K = 2, K' = 1. For the A-if setup, set 
A 1 = A, A 2 = t, Li = L and L 2 = /; for the t-if setup, 
set A 1 = t, A 2 = A, L\ — I and L 2 = L. 



Appendix D: Piecewise linear potential 

We solve the Fokker Planck equation associated with 
Eq. [7] in Laplace space, 

sp -p + d x (sgn(x)aDp - Dd x p) = 0, (Dl) 

where p — p(x,s) — p(x, t)e~ st dt is the Laplace 

transformed density, and the initial condition pa(x) — 
p(x,t = 0_) = 0. Particles are injected at t = 0; 
the boundary conditions for the total flux j(x,t) = 
sgn(x)aDp(x,t) — Dd x p(x,t) at the left boundary read 

j(x = —1, t) = 8(t) or in Laplace space, 
j(x = -l,s) = l. (D2) 

They incorporate both the reflecting boundary and the 
injection of unit probability ata; = — 1 at t = 0. 

At the right boundary, we consider cither reflecting 
(referred to as r/r) or absorbing (r/a) conditions: 



(B4) 



j(x = 1,8)- 

p(x = l,s) 



0, or 

= 0, respectively. 



(D3) 
(D4) 
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Eqs. [DTJ |D2] and [D3HD4 can be solved by using the 
ansatz p(x, s) — e 2<lx p(x, s), and joining solutions in the 
regions x < and x > 0. After straightforward but 
lengthy algebra, the solution in the r/r case can be writ- 
ten as 

p(x,s) = e -i a(1 - M) x 

qa sinh(q(l— — 2q 2 cosh(q(l — r)) — 9( — r)a 2 sinh(gr) sinh(g) 
2sa sinh 2 (q) — Asq sinhfg) cosh(g) 

(D5) 

where q = \J s/ D + a 2 /4 and 9 is the unit step function. 
For the r/a case we obtain 



p(x,s) = 

qa sinh(q(l — r)) 



|a(l-|r|), 



-6(— r) i L(cosh(q(l+r)— cosh(g(l-r))) 



2sa sinh 2 (q)-\-aDq 2 

. (D6) 

Both Green's functions have poles only on the non- 
positive real s-axis. The barrier crossing time can be ex- 
tracted by solving for the largest negative pole at — sab! 
in the r/r case, sab = ^ab + ^ba = 2&AB while in the r/a 
case, sab = &ab- If the barrier is high, we may expand 



D5 



D6 



for 



Da 2 



< 1. 



the relevant denominators in eqs. 

We find in both cases (r/r and r/a) that the waiting time 
scales exponentially with the barrier: 

2e a 



TAB 



'AB 



a 2 D 



0(l/a). 



(D7) 



In the r/a case we can also evaluate the exit flux through 
the absorbing boundary, 

Aq 2 



j(x = l,s)= jab(s) 



(D8) 



2 - (a 2 - V) cosh(2 9 ) 

The equilibration time within basin A can be estimated 
as a diffusion time for covering the thermally accessible 
range of x; alternatively, a more accurate pre- factor can 
be obtained by solving Eq. |Dl| with Eq. |D2| as above but 
replacing the W-shaped potential by a uniformly increas- 
ing ramp potential U = ax of the same slope. Evaluating 
the slowest relaxation time now gives the relaxation time 
in A, 



TA = 



a 2 D 



(D9) 



Finally, the crossing time tq is the mean first passage 
time for diffusion through region C; that is, from the 
boundary of region A at — ip, up the barrier and down 
on the other side until reaching the boundary of region 
B at xc, without returning to A. If we disregard trajec- 
tories that cross x — more than once, tq is the sum of 
the mean first passage times for the two segments, with 
negative and positive constant drift, respectively. The 
non-intuitive but well-known result is that diffusion of 
successful transition paths up the barrier takes as long 
as down the barrier (see e.gp3). In the drift-dominated 
regime, the mean first passage time is controlled by the 
drift velocity aD. One obtains 

2 



tc 



2 x 



< 



aD ~ aD 



(D10) 
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